package cn.edu.ncepu;

import java.util.ArrayList;
import java.util.List;

public class Matrix {
    private final List<List<Double>> lines;
    private static final double MAX_ERROR = 1e-15;

    // 构造函数
    public Matrix(List<List<Double>> lines) {
        // 深拷贝输入矩阵，避免外部修改影响内部状态
        this.lines = new ArrayList<>();
        for (List<Double> line : lines) {
            this.lines.add(new ArrayList<>(line));
        }
    }

    // 获取矩阵行（返回拷贝，避免外部修改）
    public List<List<Double>> getLines() {
        List<List<Double>> copy = new ArrayList<>();
        for (List<Double> line : lines) {
            copy.add(new ArrayList<>(line));
        }
        return copy;
    }

    // 矩阵转字符串
    @Override
    public String toString() {
        StringBuilder sb = new StringBuilder();
        for (List<Double> line : lines) {
            sb.append(line).append("\n");
        }
        return sb.toString();
    }

    // 矩阵列拼接（this | matrix）
    public Matrix append(Matrix matrix) {
        int rows = this.lines.size();
        int colsThis = this.lines.get(0).size();
        int colsOther = matrix.lines.get(0).size();
        int totalCols = colsThis + colsOther;

        List<List<Double>> unit = new ArrayList<>();
        for (int y = 0; y < rows; y++) {
            List<Double> row = new ArrayList<>();
            // 添加当前矩阵的列
            row.addAll(this.lines.get(y));
            // 添加拼接矩阵的列
            row.addAll(matrix.lines.get(y));
            unit.add(row);
        }
        return new Matrix(unit);
    }

    // 子矩阵截取 [rowStart, rowEnd) × [colStart, colEnd)
    public Matrix subMatrix(int rowStart, int rowEnd, int colStart, int colEnd) {
        List<List<Double>> unit = new ArrayList<>();
        for (int i = rowStart; i < rowEnd; i++) {
            List<Double> row = new ArrayList<>();
            for (int j = colStart; j < colEnd; j++) {
                row.add(this.lines.get(i).get(j));
            }
            unit.add(row);
        }
        return new Matrix(unit);
    }

    // 矩阵减法
    public Matrix subtract(Matrix other) {
        int rows = this.lines.size();
        int cols = this.lines.get(0).size();
        List<List<Double>> result = new ArrayList<>();

        for (int y = 0; y < rows; y++) {
            List<Double> row = new ArrayList<>();
            for (int x = 0; x < cols; x++) {
                double val = this.lines.get(y).get(x) - other.lines.get(y).get(x);
                row.add(val);
            }
            result.add(row);
        }
        return new Matrix(result);
    }

    // 矩阵乘法（this × other）
    public Matrix multiply(Matrix other) {
        int rowsThis = this.lines.size();
        int colsThis = this.lines.get(0).size();
        int colsOther = other.lines.get(0).size();

        List<List<Double>> result = new ArrayList<>();
        for (int i = 0; i < rowsThis; i++) {
            List<Double> row = new ArrayList<>();
            for (int j = 0; j < colsOther; j++) {
                double sum = 0.0;
                for (int k = 0; k < colsThis; k++) {
                    sum += this.lines.get(i).get(k) * other.lines.get(j).get(k);
                }
                row.add(sum);
            }
            result.add(row);
        }
        return new Matrix(result);
    }

    // LUP分解：返回 [L, U, P]
    public Matrix[] lupDecomposition() {
        int n = this.lines.size();
        // 初始化L、U、A为原矩阵的深拷贝
        List<List<Double>> l = deepCopy(this.lines);
        List<List<Double>> u = deepCopy(this.lines);
        List<List<Double>> a = deepCopy(this.lines);
        // 置换矩阵的索引数组
        int[] p = new int[n];
        for (int i = 0; i < n; i++) {
            p[i] = i;
        }

        for (int i = 0; i < n; i++) {
            // 寻找主元行（最大绝对值）
            int swapLine = i;
            for (int k = i + 1; k < n; k++) {
                if (Math.abs(a.get(k).get(i)) > Math.abs(a.get(swapLine).get(i))) {
                    swapLine = k;
                }
            }
            // 检查主元是否为0（考虑精度）
            if (Math.abs(a.get(swapLine).get(i)) <= MAX_ERROR) {
                throw new RuntimeException("矩阵无解：主元为0，矩阵奇异");
            }

            // 交换A的行
            swapRows(a, i, swapLine);
            // 交换L的行
            swapRows(l, i, swapLine);
            // 交换置换索引
            int temp = p[i];
            p[i] = p[swapLine];
            p[swapLine] = temp;

            // LU分解核心逻辑
            l.get(i).set(i, 1.0); // L的对角线为1
            u.get(i).set(i, a.get(i).get(i)); // U的对角线为A的主元

            for (int j = i + 1; j < n; j++) {
                // L的上三角部分置0
                l.get(i).set(j, 0.0);
                // 计算L的下三角元素
                l.get(j).set(i, a.get(j).get(i) / a.get(i).get(i));
                // U的下三角部分置0
                u.get(j).set(i, 0.0);
                // U的上三角元素等于A的对应元素
                u.get(i).set(j, a.get(i).get(j));

                // 计算舒尔补，更新A的子矩阵
                for (int k = i + 1; k < n; k++) {
                    double val = a.get(j).get(k) - l.get(j).get(i) * a.get(i).get(k);
                    a.get(j).set(k, round(val, 2)); // 保留2位小数，与Python一致
                }
            }
        }

        // 构建置换矩阵P
        List<List<Double>> pMatrix = new ArrayList<>();
        for (int i = 0; i < n; i++) {
            List<Double> row = new ArrayList<>();
            for (int j = 0; j < n; j++) {
                row.add(p[i] == j ? 1.0 : 0.0);
            }
            pMatrix.add(row);
        }

        return new Matrix[]{new Matrix(l), new Matrix(u), new Matrix(pMatrix)};
    }

    // 解线性方程组 Ax = b（通过LUP分解）
    public List<Double> solveByLup(List<Double> values) {
        Matrix[] lup = this.lupDecomposition();
        Matrix l = lup[0];
        Matrix u = lup[1];
        Matrix p = lup[2];

        // 构建b的列矩阵
        List<List<Double>> bMatrix = new ArrayList<>();
        for (Double val : values) {
            bMatrix.add(new ArrayList<>() {{
                add(val);
            }});
        }
        Matrix pbMatrix = p.multiply(new Matrix(bMatrix));
        // 提取Pb的数值
        List<Double> pbValues = new ArrayList<>();
        for (List<Double> row : pbMatrix.getLines()) {
            pbValues.add(row.get(0));
        }

        // 调用solve方法求解
        return solve(l, u, pbValues);
    }

    // 辅助方法：解 Ly = Pb，再解 Ux = y
    private List<Double> solve(Matrix l, Matrix u, List<Double> values) {
        int n = values.size();
        List<Double> y = new ArrayList<>(values); // 初始化y = Pb

        // 前向替换解Ly = Pb（L是单位下三角矩阵）
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                double val = values.get(j) - l.getLines().get(j).get(i) * values.get(i);
                values.set(j, val);
            }
        }

        // 后向替换解Ux = y
        List<Double> solution = new ArrayList<>();
        for (int i = 0; i < n; i++) {
            solution.add(0.0);
        }
        for (int i = n - 1; i >= 0; i--) {
            double val = round(y.get(i) / u.getLines().get(i).get(i), 2);
            solution.set(i, val);
            for (int j = 0; j < i; j++) {
                double newY = y.get(j) - u.getLines().get(j).get(i) * solution.get(i);
                y.set(j, newY);
            }
        }
        return solution;
    }

    // 辅助方法：深拷贝二维列表
    private List<List<Double>> deepCopy(List<List<Double>> original) {
        List<List<Double>> copy = new ArrayList<>();
        for (List<Double> line : original) {
            copy.add(new ArrayList<>(line));
        }
        return copy;
    }

    // 辅助方法：交换二维列表的两行
    private void swapRows(List<List<Double>> matrix, int i, int j) {
        List<Double> temp = matrix.get(i);
        matrix.set(i, matrix.get(j));
        matrix.set(j, temp);
    }

    // 辅助方法：保留指定位小数
    private double round(double value, int places) {
        if (places < 0) throw new IllegalArgumentException();
        long factor = (long) Math.pow(10, places);
        value = value * factor;
        long tmp = Math.round(value);
        return (double) tmp / factor;
    }
}